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Abstract 

The construction of multigrid operators for disordered linear lattice operators, in particular the 
fermion matrix in lattice gauge theories, by means of algebraic multigrid and block LU decompo- 
sition is discussed. In this formalism, the effective coarse-grid operator is obtained as the Schur 
complement of the original matrix. An optimal approximation to it is found by a numerical op- 
timization procedure akin to Monte Carlo renormalization, resulting in a generalized (gauge-path 
dependent) stencil that is easily evaluated for a given disorder field. Applications to preconditioning 
and relaxation methods are investigated. 
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I. INTRODUCTION 

Most of the computer time in lattice gauge theory simulations is spent today on inverting 
the Dirac matrix, which describes the dynamics of quarks. While this problem can be tackled 
efficiently by iterative solvers such as Krylov subspace methods |2|, |3|, these solvers suffer 
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from critical slowing down in the physically interesting regions of small quark masses. One 
of the limiting factors here is that the elementary operation in iterative solvers, the matrix- 
vector multiplication, affects only next neighbors and thus limits the speed of information 
propagation over the lattice (though it can be improved by ILU preconditioning 0, EJ). Using 
multigrid methods, one attempts to overcome this limitation by propagating information on 
a hierarchy of coarser and coarser grids. 

Multigrid methods || [7| for the Dirac matrix inversion were inspired by Mack's multigrid 
approach to quantum field theory |§ and have been proposed and investigated around 1990 



by groups in Boston gg, y, [L2|, |13|, 0J, Israel p|, 0, y, |J, Amsterdam pj, |0|, H |2 



and Hamburg |23|, |24|, |25|, |26], |27], |28| . Though it was shown that multigrid methods would, 
in principle, greatly reduce or eliminate critical slowing down, they have not been able 
to improve the performance of the Dirac matrix inversion in actual applications. These 
classical multigrid methods are based on a geometrical blocking of the lattice that leads 
to an effective coarse-grid formulation of the original matrix by a Galerkin approach. The 
choice of the blocking and interpolation kernels is crucial for the performance of the multigrid 
algorithm; they must be chosen such that the long-range and short-range dynamics decouple 
as much as possible. An optimal choice of the blocking kernels thus depends on the gauge 
field and must in principle be recalculated for each gauge configuration as the solution of a 
minimization problem. As this is not feasible, various approximations have been introduced 
in the literature. 

Since then, the algebraic multigrid method [|(], [31], [32|, [33|, [34[] has emerged as a proven 
new method in numerical mathematics which does not rely on a geometrical decomposition 
of the lattice, but solely on the matrix itself. It starts out from a thinned lattice that contains 
a certain subset of the original lattice sites retaining their original values and constructs an 
effective operator on this lattice by means of an block LU decomposition. Thus, instead of 
geometric proximity, the actual matrix elements determine in the coarsening procedure. For 
lattice gauge theory, this would automatically take the dynamics of the gauge field, which 
appears as a phase factor in the matrix elements, into account. (A similar reasoning was 



already cited by Edwards, Goodman, and Sokal in one of the first papers on multigrid 
methods for disordered systems, namely on random resistor networks, though it was later 
assumed that algebraic multigrid would be much more costly than the geometric variant 



TTfl.) Recently, algebraic multigrid was applied to lattice gauge theory by Medeke [36 
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and conversely research about renormalization-group improved and perfect operators in 
lattice gauge theory has been used to obtain efficient coarse-graining schemes for partial 
and stochastic differential equations |J7 . 

One advantage of algebraic multigrid over geometric methods is that it directly gives 
a prescription for calculating the coarse-lattice operator and the interpolation kernel by a 
submatrix inversion problem. For this purpose, it is instructive to consider the problem 
of matrix inversion not in terms of the spectral decomposition, but in terms of the von 



Neumann series [38, 39]. Combined with the block LU decomposition, this results in a simple 
computational picture of the method as a noninteracting random walk on a hierarchically 
decomposed lattice in a disordered background. The resulting algebra of paths on this 
lattice is the basis for the computational construction of multigrid operators below, where 
the expansion coefficients are calculated numerically in a process similar to Monte Carlo 
renormalization. In that way, physical (or rather, computational) information about the 
fine-lattice operator is harvested to approximate the coarse-lattice operator as closely as 
desired. 

The lattice fermion matrix is an example of the more general problem of disordered 
matrices operators which has many applications in different fields of computational physics. 
Some recent instances in the literature are e.g. the transport properties of light particles 
in solid-state physics jll], [01], or transport in porous media, where multigrid methods have 
already been employed |4*2| , ?3|. The numerical method investigated here is sufficiently 
general to be applicable to such problems. 

Another factor that might turn out in favor of multigrid algorithms comes from recent 
advances in computer technology. Supercomputer architecture has long been dominated by 
vector-oriented machines that can quickly execute relatively simple and uniform elemen- 
tary operations on large vectors. But because processor speed advances much faster than 
memory access speed, the performance of such operations is now severely limited by memory 
bandwidth. In computer architecture, this is overcome by using cache-oriented architectures 



such as in modern microprocessors, from which e.g. cluster computers |44| are built. These 
architectures are efficient for algorithms that have a high balance, i.e. number of operations 
per memory access, and make a frequent reuse of memory locations [[45] , SSL but can tolerate 
much more complexity and irregularity. Multigrid algorithms usually result in denser, but 
smaller matrices and might thus be favored by these architectures. 
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Finally, problems other than matrix inversion have recently come into focus for the Dirac 
operator. One is the calculation of its low-lying eigenvalues, which was one of the first appli- 
cations of multigrid [|27] and is today used for investigating hadron structure j|7| and chiral 
properties [58, 55]. Another one is the development of chirally-improved Dirac operators 
such as the Neuberger operator |5t], 5TJ which requires the calculation of the inverse square 
root of a matrix. Algorithms for this problem are still under development, but as they 
also employ iterative solvers, it can be hoped that multigrid operators can aid in improving 
performance there. Multigrid methods have also been cited as a possible improvement to 



the domain-wall formulation of lattice gauge theory on a five-dimensional lattice |p2 



The layout of the paper is as follows: In section 2, we introduce the lattice Klein-Gordon 
and Wilson-Dirac operator and their interpretation as diffusion on a disordered lattice, and 
in section 3, the general algebraic multigrid procedure. The numerical procedure to construct 
multigrid operators is discussed in section 4, and some applications to preconditioning and 
relaxation algorithms are discussed in section 5. 



II. MODEL 



The model considered here is the Euclidean lattice Klein-Gordon and Dirac operators with 
a (7(1) gauge field, discretized on a square lattice of lattice spacing a = 1. The formalism 
will also apply to other gauge groups and can be generalized to other operators and other 
types of disorder. Numerical calculations are performed in two dimensions and on relatively 
small lattices (16 x 16), so that the complete spectrum of the operators can be analyzed. 



A. Klein-Gordon operator 

For simplicity, we start with the Klein-Gordon operator which, without a gauge field, is 
conventionally written as 

- A(x,y) +m 2 5 X) y (1) 

with the next-neighbor lattice Laplacian A(x, y) and the (bare) mass m. Its eigenmodes 
are discrete Fourier modes, and its Green's functions decay exponentially as e _m ' x ~ y ' with 
the inverse correlation length m. Gauge disorder is introduced as a U(l) gauge phase U^(x) 
associated with the links (x, fi) of the lattice. The disordered Klein-Gordon operator can be 
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rewritten, up to a normalizing factor, as 



M(x,y) = 5 X)V - KQ{x,y) , 

±D 

Q{x,y) = $3<We M E^(aO • ( 2 ) 
11=1 

Here x, y are sites on the lattice, \i labels positive and negative directions ±1, . . . ,±D, 
is the unit vector in direction //, and we use the convention = — e M and U-^x) = 
U^(x — e M ). The completely hopping parameter k are borrowed from fermion terminology 
and related to the bare mass m in (|l|) by 

1 



K 



(3) 



m 2 + 2D 

It determines the spectrum of M and thus the physical properties of the theory. For the free 
theory, i.e. U^(x) = 1, the eigenvalues of Q lie between —2D (checkerboard configuration) 
and 2D (completely constant configuration). Consequently, the eigenvalues A n of M lie in a 
band of width 2Dk around 1: 

1 - 2Dk < X n < 1 + 2Dk . (4) 
For values of k below a critical value 1/2D, 

K< ^ = Kcrit ( 5 ) 

the spectrum of M is strictly positive and M invertible. The situation for a theory in a 
gauge field will be discussed in the next section. 

Physically, the interesting quantity in M is the location of the lowest eigenvalue that 
determines the mass gap of the theory and the exponential decay of its Green's functions, 
i.e. the correlation functions or propagators. As k approaches k ci h, the mass gap becomes 
smaller and the correlation length increases, until at t\i — ^crit 5 ^ Z6rO 61 genvalue appears 
and the correlation functions, after proper subtraction, decays logarithmically or rationally, 
depending on the dimension D. 

The performance of iterative solvers for the inverse of M is linked to the condition number, 
i.e. the ratio of the largest to the smallest eigenvalue, that diverg Unfortunately, 
this is exactly the region that is of interest in lattice gauge theories, where quarks are much 
lighter than the typical scales of the theory. 
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B. Diffusion representation of M 



Conventionally, eq. (0) is interpreted in lattice field theory in terms of its spectrum that 
determines the dispersion relation and thus the mass of the physical particle. However, this 
description is less useful in the presence of a gauge field, as there is no simple expression for 
the eigenmodes and eigenvalues of M. Alternatively, M can also be interpreted as describing 
the diffusion of noninteracting random walkers, i.e. Brownian motion, on the lattice. This 
can most easily be seen by using the von Neumann series to calculate M _1 : 

(1-kQY 1 = 1 + kQ + (kQ) 2 + (kQ) 3 + ... 

oo ±D ±D 

(1-kQ)-£ = S y,x+e n +-+e l > n K n V(U;x,n 1 ,...,n n ) . (6) 

n=0 Mi 

The sum here runs over all connected paths leading from x to y, and each path carries a 
gauge phase assembled from the links it traverses of the lattice: 

Pi*i,...,i* m (U;x) = V(U;x,x + e fll ,...,x + e lll + ... + e fln ) (7) 
= U IM1 (x)U IM2 (x + e IMl )...U tln (x + e tll + --- + e IMn _ 1 ) . (8) 

The diffusion interpretation is based of the fact that the sum in eq. can be generated by 
a random walk starting at x and making random moves along the links of the lattice until 
reaching site y. Since the probability for each path of length n is (2D)~ n , eq. (||) can be 
rewritten as an ensemble average over random walkers: 

where the average is over all random paths of length n, and each random walker carries a 
phase factor. 

In the free case, and for k = 2D, eq. (^) is just the time-integrated probability, or average 
time spent, at site y for a random walker starting out at x. Other values of k correspond 
to spontaneous creation or absorption of walkers. In the presence of gauge disorder, each 
walker additionally carries a phase factor V(U) picked up from the links it has travesed. 
Their contributions can therefore interfere destructively, leading to a stronger decay of the 
Green's function than in the free case: 

V-kQ{U))-1 < (1 — kQ(U = O))^ 1 . (10) 



In particular, this means that for those k for which the free correlation function exists, 
i.e. for k < 1/2D, the disordered correlation function also converges, and the critical hopping 
parameter K crit of the disordered theory therefore must be equal or larger than 1/2D. From 
the discussion in the preceding section, this reflects an upward shift of the lowest eigenvalue 
of M that is interpreted as a dynamically generated mass. 

The value of the dynamical mass depends on the correlation characteristics of the gauge 
field. If the gauge field is strongly correlated, similar paths will have similar phase factors, 
the destructive interference will be less severe and the dynamical mass smaller than for more 
disordered fields. In this way, the diffusion representation provides a qualitative picture of 
the spectrum of M in a disordered field. 



C. Wilson-Dirac operator 

The Dirac operator describes the fermionic quark fields in lattice gauge theory simu- 
lations. These fields possess an internal spin degree of freedom and are described by a 
[d/2 J -component spinor field. The Dirac equation is a first-order differential equation, and 
its naive discretization on the lattice suffers from fermion doubling, which is removed in the 
Wilson-Dirac discretization: 

M D (x,y) = S x , y - K,Q D (x,y) , 

±D 

Q D (x,y) = Yt^x+eAi-^U^x) . (11) 
where in two dimensions the Dirac matrices 7^ are given by 

Together with the unit matrix and the matrix 75 = 27172, they form a basis for the linear 
operators in spinor space. The Wilson-Dirac operator is a non-Hermitian operator and thus 
has in general complex eigenvalues and different left and right eigenvectors, but a Hermitian 
Wilson-Dirac operator can by defined by multiplying it with 75: 

M HD = l5 M D . (13) 

For the matter of inversion, this operator is completely equivalent to the original Wilson- 
Dirac operator. Due to the two-component spinors it acts upon, it has twice as many degrees 



of freedom as the Klein-Gordon operator. These additional degrees of freedom show up in 
the spectrum of Mhd as a band of negative eigenvalues mirroring the positive band seen 
in the Klein-Gordon operator. The argument used above about the relation of k to the 
low-lying eigenvectors and the mass gap remain valid for the Dirac operator. 



III. ALGEBRAIC MULTIGRID 



A. Schur complement 

We consider the solution of a linear system 

J2M(x,y)f(x) = a(y) (14) 
y 

with the right-hand side a given and / unknown. In the algebraic multigrid approach, the 
sites of the lattice £ are decomposed into a coarse lattice C\ and a fine complement lattice 
C 2 = £\£i. This decomposition can in principle be performed without reference to the 
geometric locations of the sites. In particular, the coarse-lattice sites are not block averages, 
but retain the same values as they had on the fine lattice. Eq. ( |14|) can then be rewritten as 

/ 



M 2 i M 22 



h 




(15) 



where the column vectors are reordered so that fi and cti are defined on the coarse lattice, 
and f2i 0,2 on the fine lattice. Applying Gaussian elimination, the fine-lattice field f2 is 
eliminated completely from the first row: 

(M11- M l2 M 22 l M 2l ) fi = (ai - M X2 M 22 l a 2 ) (16) 
M 22 f 2 = (a 2 -M 21 f 1 ) (17) 

The matrix operator in front of f\ plays the role of an effective operator on the coarse lattice. 
It is called the Schur complement Sm of the matrix M with respect to the decomposition 

u 

S AI = Mn - M,.\/, 2 ; M M . (18) 

The simplest example of a Schur complement is the odd-even decomposition of a square 
matrix with next-neighbor interactions, where M 22 is diagonal. 
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Since the Schur complement contains the inverse of the fine-grid submatrix M22, algebraic 
multigrid must strive to find decompositions that make M22 as dominated by the diagonal as 
possible so that a good approximation to its inverse can be found. The method thus provides 
an immediate heuristic for choosing a multigrid decomposition based only on the content of 
the matrix M which is advantageous e.g. for irregular discretizations. In the case of lattice 
gauge theories, the coefficients fluctuate stochastically on the gauge group, but have the 
same magnitude. We thus choose a regular decomposition where C\ = (2Z) D , i.e. the coarse 
lattice is the lattice of points with all-even coordinates (cf. Fig. |l|; for a different choice see 
j3Gfl). Since C\ is not connected with respect to M, M\\ is diagonal, M i2 links points of C\ 
to the remaining points, and M 2 2 forms a connected lattice on the fine grid. 
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FIG. 1: Two-grid decomposition of a two-dimensional lattice. Filled points are lattice sites of the 
coarse grid (£1), open points sites of the fine grid (£2). Dashed lines are links that enter in Myi 
and M21, while solid lines are those in M22- The arrows show three possible contributions of order 
2, 4, and 6 to a general coarse-grid operator of the form (Hjj). 
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B. Block LU factorization 



Eq. (|16D and flT7|) can be expressed in matrix form as an block LU decomposition of M. 
Its general form is 

M (i-Q + ) (so\( i o 

^0 i ){o T) {-Q 1 

The advantage of the representation ( |IT| ) is that its inverse can be readily written as 



(19) 




M _ 1= „ . _, " (20) 



requiring only the inversion of the smaller matrices S and Q. In multigrid terminology, 
applying the factorization fl20p to a vector amounts to first applying an restriction operator 
Q + that moves information from the fine grid to the coarse grid, then solving the effective 
coarse-grid operator S along with the residual fine-grid operator T, and finally using the 
interpolation operator Q to move information back from the coarse grid to the full grid. 
Inserting (|l9|) into flTljD leads to the following relations that define S, Q, and T. 

Sfi = ai -Q + a 2 (21) 
T{-Qh + h) = a 2 . (22) 

Comparing with eq. (|T6|) and ([17]) gives the solution 

S = S M = M n - M 12 M^ 2 1 M 21 , 
Q = -M^M 2l , 

T = M 22 . (23) 

The Schur complement is the exact effective coarse-grid operator associated with a given 
decomposition of the lattice. Unlike the geometrical approach, this decomposition is per- 
formed simply by thinning the matrix and thus without a block averaging that might in- 
troduce artifacts in the presence of a gauge field, and the resulting coarsened operator is 
defined solely in terms of the original matrix and thus the original disorder field. In the 
following sections, we will show that this procedure can be interpreted in terms of a nonin- 
teracting random walk leading to an explicit representation of the effective matrix in terms 
of a generalized path-dependent stencil. 
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C. Renormalization group and projective multigrid 



The Schur complement has a close relation to the renormalization group approach used 
in statistical mechanics, where the Klein-Gordon and Dirac fields are described by a path 
integral with a Gaussian weight. In this approach, the partition sum 

df exp [-!(/, M/) + i(a,/) + c.c." 
= const • e^ M ~ la) 

(the integration is over all lattice field configurations /) generates the inverse of M by taking 
the derivative with respect to the source field a(x): 

1 d 2 



Z da*(x) da(y) 



\M-\ Xl y) . 



(24) 



a=0 



In renormalizat ion-group inspired multigrid approach, such as proposed by Mack ||, 
Brower et. al. |12j or Vyas [[TJJ, one introduces new variables in the path integral by means 
of coarsening. In our terminology, a coarsened lattice field Fi defined on on £1 is given by 



Fx = h + Cf 2 



(25) 



where the coarsening matrix C performs an averaging over neighboring sites in £2- (A more 
general coarsening would also include an averaging over neighboring sites from Ci). The 
coarsened field F\ is interpolated back to the fine lattice using an interpolation kernel A and 
subtracted from f 2 to form the residual fine-grid field: 



C2 = h ~ AF, 



(26) 



In this way, it is hoped that some of the dynamics of the fine lattice £2 is moved to the 
coarse lattice. The whole transformation reads 

(27) 




V 



1 

-A 1 




As its determinant is unity, F\ and (2 can be used as new integration variables in the path 
integral. The quadratic form in the exponential is then 

F 1 \ + 1 1 A + \ I 1 
C2 ) \0 1 ) \A1 

-c+ 1 01 



(/, Mf) 



M 




(28) 
(29) 
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and the coupling to the external fields 



I \ + 

' ai + A + a 2 - A + C + a 1 




(30) 



, d 2 — C CLi 

Algebraic multigrid does not use coarsening; the coarse lattice field is obtained by thinning 
out the original lattice and thus F\ = fx, C = and M = M. Then the block LU 
decomposition (19) tells us that we can make the quadratic form ( p8|) block diagonal by 



choosing A = — M 22 1 M2i, in which case the partition sum factorizes into a coarse-grid and 
a fine-grid integral: 

r 1 1 /* 1 1 

2 = / e -2(/i' S 'M/i)+ 2 -(ai+A+a2,/i)+c.c. _ / {fl,M22h)+2 (o2,/2)+c.c. ^g-Q 

The coarse scale dynamics is obtained by taking derivatives of Z with respect to a\ at 
a2 = and is therefore completely contained in the first integral, which is governed by Sm, 
the Schur complement. Thus the algebraic multigrid choice of the interpolation kernel A 
can be characterized as that choice that leads to a complete decoupling of fine and coarse 
degrees of freedom, which makes integrating out the fine-grid degrees of freedom £ 2 trivial. 
Note that the Schur complement can be used to obtain effective fermions such as used in 
the blocked fermion approach to full QCD |53 |. 



If a nonvanishing averaging kernel C is used, one can achieve complete decoupling in the 
quadratic form by applying the Schur complement to the matrix M, but now a\ also couples 
to (2 in the source term (P0|), and the integral over (2 results in a nontrivial contribution 
to the partition sum. This contribution is governed by M 2 ~ 2 1 , the propagator of the residual 
modes ( 2 on the fine lattice, which is given by 

M 22 = M 22 - M 21 C - C + M 12 + C + M n C . (32) 

One therefore strives to choose an averaging kernel C such that M22 1 is as local as possible, 
indicating that most of the (long-range) dynamics has been moved to the coarse-lattice field 
F x . 

The decoupling of coarse and fine modes was at the heart of Mack's original multigrid 
proposal H , which was inspired by constructive quantum field theory. It was already noticed 
by Brower, Rebbi, and Vicari |TU| that minimizing the coupling between coarse and fine 



modes amounts in spirit to algebraic multigrid, and using a similar reasoning, the form of 
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the Schur complement is mentioned in eq. (12) of Jl2| by Brower, Edwards, and Rebbi. In 
practice, in the projective multigrid method the coarsening C was chosen to localize M22 
as much as possible, and the interpolation A to maximize the decoupling between coarse 
and fine modes. An explicit method for calculating A for a given C is specified as idealized 
multigrid algorithm by Kalkreuter |27J , and it can be verified that for a "lazy" coarsening 
kernel C = 0, the resulting interpolation kernel has the form (^). The actual calculation 
for nontrivial C is numerically infeasible as it requires the solution of a nonlinear equation 
at each coarse lattice site, so instead the ground-state of a block-truncated matrix was used 
for the kernels, leading to the method known as ground-state projection multigrid. The 
effective coarse-grid matrix was then approximated by a Galerkin choice using the rows of 
A as a basis. 

Using the algebraic multigrid choice of C — 0, we cannot optimize the decay of M^ 1 . 
However, as there is no induced coupling of ai to the fine lattice field ( 2 , the Schur com- 
plement provides an explicit representation of the complete coarse-lattice dynamics. In the 
following, we will discuss the interpretation of this quantity based on the diffusion represen- 
tation. 



D. Diffusion on a multigrid 



In sec. [II B| , the diffusion representation of the inverse of a matrix M was introduced. 
The matrix elements of the Schur complement Sm define a similar effective diffusion process 
involving only coarse lattice sites with transition probabilities and gauge phases determined 
by the elements of Sm- As Sm is identical with the upper left submatrix of M _1 , 

s m( x ,v) = M ~ l (. x >y) tor x,y e d , (33) 

this diffusion process is equivalent to the original diffusion process when origin and destina- 
tion of the random walk are restricted to the coarse lattice. Using (|18D, the matrix elements 
of Sm can itself be written in a diffusion representation: 

S M (x,y) = 8 x , y - Yl K n V(U\x,xi, . . . ,x n ,y) . (34) 

n 

(x, X\ , . . . , x n , y) 

X\, . . . , X n E L-2 
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Since the sum comes from expanding M^ 1 , the inner parts (xi, . . . ,x n ) of the paths run 
exclusively on the fine lattice £ 2 ; they are connected to the coarse lattice C\ at points x and 
y through the matrix elements of M 12 and M 21 (some of these paths are shown in Fig. [[]). 
The matrix elements of the Schur complement can therefore be interpreted as resulting 
from random walks on the fine complement lattice £2, and when the Schur complement 
is inverted, it defines a random walk on the coarse lattice C\, in each of whose steps all 
possible paths on the fine lattice £2 are effectively taken into account. In this way, the 
Schur complement allows to "integrate out" the fine degrees of freedom in the diffusion 
process much as renormalization group transformations do. 
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FIG. 2: Two-grid decomposition of a odd-even decomposed lattice. As in Fig. [I], filled circles 
mark sites of the coarse grid, and open circles those of the fine grid. Full lines show links from 
M22, dashed lines those from Mix, and dotted lines those from Myi. M\\ and M22 here form two 
similar grids, with M\% connecting the two along the diagonals. 

In the two-dimensional case, the situation can be made more symmetrical by performing 
and even-odd-decomposition before the multigrid decomposition. In an even-odd decompo- 
sition, the matrix is replaced by the effective matrix on the even sites of the lattice, 

S ee = 1 - M eo M oe , (35) 
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where M eo and M oe are the matrix elements of M between even and odd sites, resp. In the 
resulting lattice topology, each site is linked both to its four next neighbors in the straight 
direction as well to the four next neighbors in the diagonal (cf. Fig. |2|). Applying now the 
multigrid decomposition, the lattice decomposes into two similar square lattices. Mn and 
M 22 are represented by the straight links on each lattice, while M 12 and M 21 are carried by 
the diagonals linking the two lattices. The effective hopping parameter on each lattice is 

2 

K e g = 2 and thus smaller than on the original lattice, and consequently Mn and M 22 
are less critical than M. In the diffusion interpretation, even if probability was conserved 
on the original lattice, it is not conserved on the two sublattices individually, as the random 
walkers can move from one lattice to the other along the diagonal links represented by M 12 
and M 21 . The criticality of the system is only regained when Mn and M 22 are combined in 
the Schur complement. 

The propagator on the fine lattice M^ 1 will therefore be short-ranged and can be incor- 
porated approximately by considering only a limited number of hops on the fine lattice in 
the expansion (0). In the following section, we exploit this to construct an approximate 
Schur complement numerically. 



IV. CONSTRUCTION OF THE MULTIGRID OPERATOR 



The central problem in deriving a multigrid algorithm for disordered systems is the con- 
struction of a suitable approximation to the Schur complement Sm an d the block LU de- 
composition (|19|). 



A. Numerically optimized Schur complement 

Our method is based on numerically approximating the Schur complement as a linear 
combination of suitable basis operators. The full Schur complement S usually cannot be 
computed exactly, but it can be characterized numerically by a sufficiently large set of pairs 
{/W, oW; n — 1, . . . , N} which satisfy 

M(U (n) )f {n) = a (n) , (36) 

where it is explicitly indicated that M depends on different realizations U^ n > of the disorder 
field. Using the block LU decomposition fllTf ), this implies that S, Q, and T satisfy the 
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relations (21), (p2|). In particular, if the pairs are chosen such that a 2 = 0, their coarse- 
lattice components {fx n \o>x } characterize the Schur complement by the relation: 

Sfi n) ^aP . (37) 

Similarly, the interpolation operator is characterized by the pairs {fi n \ f^} using 

Qf[ n) = A n) ■ (38) 

There is no need to approximate T = M 22 , as it does not contain the inverse of M 22 . An 
approximation S to S, and similarly Q to Q, can then be characterized with minimum bias 
by the mean error norm of (|37]): 

N 2 

6 2 = J2\S(U^)f[ n) -a^ (39) 

n 

This quantity is an approximation to an operator norm \SS — 1| 2 averaged over gauge field 
configurations. The choice of pairs (|36|) introduces a weight function in this norm and thus 
determines which functions are well approximated by S. Eq. (|39|) can be rewritten as 

N N 

S 2 = \(S(U^)S - l) 4 n) | = - S(U (n ^) ft ] \ (40) 

n n 

If fx is chosen randomly on the whole function space, 5 2 thus approximates \S — S\ 2 , while 
if ax is chosen randomly, it approximates 

\SS- 1 -1\ 2 = \(S-S)S~ 1 \ 2 (41) 

The factor 5* _1 increases the weight of functions that have eigenvalues closer to zero, and 
these are exactly the long-ranged functions we are interested in. Choosing the right-hand 
sides ax also allows us to enforce a 2 = and thus remove Q from the (pT|) . In practice, we 
choose a to be a delta function on a randomly chosen coarse- lattice point, and determine 
/ by solving the linear system (|36D . In other words, the f^ are Green's functions, or 
propagators in lattice gauge theory language, and we look for an approximate operator S 
that well reproduces the action of S on these functions. 

If we choose a general linear combination of suitably chosen basis operators S^(U) re- 
specting the symmetries of the problem as approximate operators, 

S(U)=Y,<* i SU(U) , (42) 
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the approximation error S 2 is a quadratic form in the coefficients af. 

5 2 = E a t A H a i ~ E a * B i ~ c - c - + c ( 43 ) 

ij i 
N 

= E (S { HU (n) )f[ n \S^\U^)f[ n) ) c (44) 

n 1 
N 

Bi = E {S {l) (U {n) )fi n \a i ? ) ) c (45) 

n 1 
N 

C = E(<4 n U n) ) £l (46) 



where (ai,bi)d denotes that scalar product on the coarse grid L\. Thus the minimum of 
the approximation error and thus the optimal choice of coefficients is found simply by: 

a, = E(^%£, (47) 

j 

The procedure for finding the coefficients of the numerical approximation to the Schur 
complement is therefore as follows: 

1. Generate right-hand sides by choosing random delta functions on the coarse grid. 

2. Find the corresponding by solving (|5B|) using the original fine-grid operator M. 



3. Project and to the coarse lattice. 

4. Calculate the matrix elements ffPf) and (1151) and find a,- from 



A similar procedure is used to approximate the interpolation operator Q. This method is 
similar to Monte Carlo renormalization group as it calculates the coefficients of the effec- 
tive coarse-grid operator by applying a block-spin transformation (in our case, simply a 
projection) to an ensemble of configurations on the fine grid. Note that the procedure is 
completely general and can be used to approximate any operator, e.g. the Neuberger or a 
renormalizat ion-group improved operator, if its action on some sample fields is known. 



B. Choice of the operator basis 

The approximate coarse-grid operator S(U) is a general operator- valued function of the 
disorder field U, but the symmetries of the problem greatly restrict the possible basis op- 
erators S^(U) from which it can be constructed. To maintain gauge covariance, each 
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contribution to the matrix element M(x, y) can be factored into a gauge-field independent 
part and a product of the gauge field links along a connected path from x to y, leading to 
the following representation of M: 

M(x,y) = J2 E rn(P)S y>x+en+ ... +e ^ n V» 1 ,..., IM (U-,x) (48) 

where the sum is over connected paths from x to y, and the reduced matrix element m(P) 
does not depend on the gauge field U. Basis operators can thus be associated with the 
possible relative paths Pi = {//i, . . . , //„}: 

S^(U)(x 1 y) = T^^^.. + ^V^, m ^(U;x) , (49) 

where is independent of x, y, and U, and serves to contain a possible internal structure 
of the field, e.g. the Dirac 7 matrices. 

Other symmetries like translation, rotation, and reflection invariance as well as hermiticity 
and charge conjugation generate further restrictions on the coefficients 014. In the calcula- 
tions, our computer code automatically evaluates discrete rotation and reflection symmetry 
for all paths and constructs basis operators that are invariant under these transformations. 
This reduces significantly the number of degrees of freedom in fl43l) . Hermiticity places con- 
straints on the real and imaginary parts of the coefficients; it is not enforced and can be used 
as a check of the accuracy of the constructed operator. The general problem of construct- 
ing generalized Dirac operators was recently brought up in the context of chirally-enhanced 
operators in 0, |5^] and discussed in general in |56 . 



The representation (|48"D has the same form as the diffusion representation fl34|), with the 
additional restriction that only paths moving exclusively on the fine lattice contribute there. 
The optimization process ( |S5| ) thus results in assigning weights to the different paths in 
the diffusion representation that reflect how much a path contributes, on average, to the 
propagator. If a complete sum over all lengths n is taken, the correct weights are of course 
K~ n , as in the von Neumann series. However, if the sum is truncated, the weights are modified 
in the optimization process and the contributions of the individual paths renormalized to 
take the obmitted paths into account. 

The basis (|9|) can be compared to the situation without a gauge field. In the terminology 
of numerical analysis, a translation-invariant operator is represented by a stencil S(a) using 
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the expression corresponding to (|48|) 

M(x,y)=J2S(a)5 y , x+a (50) 

a 

where a runs over the possible offsets between lattice points. The basis operators in this 
case can therefore be labeled by the offset a, and the equation corresponding to ( pEUp reads 

S®(x,y)=6 y , x+ai . (51) 

Comparing the two cases, one can see that in the disordered case the relative paths take 
over the role of the offsets, and a generalized stencil is given by assigning each relative path 
a certain weight. The optimization procedure corresponds to calculating such a generalized 
stencil, much as an ordinary stencil is constructed for differential equations. 

Once the disordered stencil is known, the coarse-grid operator 5 for any gauge field can be 
constructed simply by evaluating all paths in the set {S^} and combining them according 
to the weights «j. Note that it is not necessary to perform the construction of the on the 
actual set of gauge fields, but on a representative ensemble only, so the optimization can be 
performed in advance of the simulation. 

The resulting coarse-grid operator contains in particular also contributions beyond direct 
neighbors on the coarse grid, similar to improved and perfect operators that lose locality but 
gain a more faithful representation of the continuum operator, and the amount of nonlocality 
can be tuned by selecting the paths in the basis set. On the other hand, if the basis set 
is chosen to reproduce the original topology on the coarse lattice, the coarse-grid operator 
defines a new effective disorder field on the coarse-grid. However, this disorder field is 
obtained by a sum and will therefore lie outside the gauge group, though it can be projected 
back to the gauge group, as is customarily done in many renormalization schemes. In this 
way, the Schur complement can be used to define a block-spin transformation for the gauge 
field. 



C. Construction in the diagonal subspace 



The basis fl49|) is the most general set of operators that can make up a gauge-covariant 
operator. Its size grows exponentially with the order of the path set and thus makes both the 
evaluation of the coefficients (44) and (|45|) and the actual construction of the approximate 
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Schur complement very time consuming. A smaller and thus more tractable basis can be 
found by considering the operators that appear in the von Neumann series (|6]) which are: 

S<® = M n , 

= m 12 M%M 21 . (52) 

This basis generates exactly the same paths as the diffusion representation ([49]) restricted to 
the fine lattice, but uses coefficients that only depend on the length of the paths. Obviously, 
if an infinite set of basis operators is used, the optimal coefficients are ±1 as given by the 
von Neumann series. However, in a finite set, the optimization procedure amounts to finding 
a polynomial approximation to the inverse function. 

To see this, let T be an approximation to M 22 , and write the approximate Schur com- 
plement as 

S = M u - M 12 TM 21 (53) 
Inserting this into (|39|) gives after some algebra 

S/i-ai| 2 =(Tf,M 21 M 12 Tf) 2 -2^(b,Tf) 2 + (b,b) 1 (54) 
with the notation (...,.. .) 2 for the scalar product over C 2 and the definitions 

/ = M 2l h , (55) 
b = Mnfx - ax , 
b = M 21 b 

= M 2l M 12 M 22 l M 2l h , (56) 

where we used fl3~7D and (|18[) in the last line. If T is a polynomial in M 22 , it is diagonal in 
the eigenbasis of M 22 and can be written as 

M 22 = Y. X kV k ®v+ (57) 

k 

T = Y.T k v k ®v+ (58) 

k 

where X k , v k are eigenvalues and -vectors of M 22 , and T k the diagonal elements of T in 
the eigenrepresentation. Inserting this into (0) gives the following quadratic form for the 
diagonal components T k of the operator T: 

Sfi ~ aif = £T* (ftfi) (MflMu^T, (59) 

kl 

-2UJ2T k (b*Jk)+ const (60) 
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where fk = (ffc, f), H — ( v k, b) are the representations of / and b in the eigenspace of M 2 2, 
and similarly for (M 2 iM 12 )h. Eq. (|56|) relates these quantities by 

h = Y,( M 2iMi 2 )kiK 1 fi (61) 

I 

which gives the following expression for (f)9|): 

Sfi - a, 2 = £ (n - i)* (A/0 ( M 2iMn) M (T, - i) + const 



2 



+ const . (62) 



Minimizing this quantity thus amounts to finding a best approximation TJ. to A^ 1 in a 
weighted square norm specified by Muf = M12M21/1. 
Choosing for T a polynomial form 

AT AT 

T=^a n+1 M| 2 " , T k = J] a n+1 A^" (63) 

n=0 n=0 

in accordance with the expansion (|52|), the approximate Schur complement then takes the 
form 

N 

S = M U -J2 a k+1 M 12 M%M 21 (64) 



fc=0 

Fig. ^| shows a numerical evaluation of such an approximation using the Dirac operator on 
a sample U(l) configuration. 

An operator of the form QS^D can be computed relatively easily: The coefficient S(x, y) 
between two sites is given by summing the phase factors over all paths from x to y on £2 with 
weights depending only on the length of the paths and given by the coefficients a k . This is a 
simplification over the full path set (|49|) , where each path can be assigned its own weight. If 
only a finite number iV of coefficients is used, the operator is local, but contains more than 
next-neighbor interactions. Note that in a typical iteration scheme the calculation of the 
matrix elements of S has to be done only once, and the complexity of the iteration depends 
only on the sparsity of the resulting S. 

D. Numerical results 

We have numerically evaluated the optimal operator both for the full path set of eq. (f49|) 
and the subset defined in eq. (|52|). The calculations were performed for the Dirac equation 
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FIG. 3: Numerical evaluation of the optimal polynomial approximation to . The curves 
show the polynomial approximations of order 1, 3, 5, and 7 to the inverse function. The actual 
eigenvalues are marked with crosses. 

on a 16 x 16 lattice in two dimensions in a U(l) gauge field generated at (3 = 3.0, which 
results in relatively smooth configurations and long-ranged correlation functions. Fig. ^ and 
[5] show sample Green's functions of the original Dirac matrix and the approximate Schur 
complement. The disorder from the gauge field shows up clearly in the Green's function, 
but the various approximations to the Schur complement reproduce the irregularities and 
the overall form of the function quite well. 

This is quantified in fig. 0, which shows the average relative error 



5 2 

(65) 



evaluated in an ensemble of ten sample U(l) gauge configurations generated at (3 = 3.0 on a 
16 x 16-lattice. For each configuration, five different Green's function at k — 0.265 on each 
configuration were calculated and the coefficients (f!4 ) and ( [45"D that characterize the error 



ellipsoid were evaluated. For the diagonal optimization, the basis (|52|) with maximum order 
i from 1 to 6, corresponding to a maximum path length from 2 to 12, was used. For the 
full optimization, a complete path basis of paths up to length 2i constrained to the discrete 
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FIG. 4: Absolute value of a sample fermionic Green's function on an 16 x 16 lattice in a U(l) 
background field for the original Dirac operator (left) and the approximate coarse-grid operator 
obtained from the Schur complement. The source is located at the origin; since the lattice is 
periodic, it shows up at the four corners of the plot. 
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FIG. 5: Cut through the example fermionic Green's function (black squares) compared to Green's 
functions calculated from an approximate Schur complement and interpolation operator of order 
1, 2, and 4. 



24 



Relative error vs. operator order 
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FIG. 6: The relative approximation error as a function of the order of the approximation, calculated 
in an ensemble of 10 independent 16 x 16 U(l) gauge configurations at (3 = 3.0 and k = 0.265. 
The full line gives the result from the optimization in the diagonal subspace, the dashed line from 
the full optimization. The dotted line is the relative error of the truncated von Neumann series of 
the same order. 
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FIG. 7: The relative inversion error as a funcjipn of the approximation order for both diagonal 
and full optimization. For comparison, the original relative error from Fig. |6| is also shown. 



rotation and reflection symmetry was taken. For comparison, the relative error resulting 
from the von Neumann series is also shown. The errors show a very clean logarithmic 
dependence on the order. As the number of paths in the operator increases exponentially 
with order, this corresponds to a polynomial dependence of the error on the number of paths. 

The full path set gives only a slightly better approximation than the diagonal approxima- 
tion, though it is significantly more difficult to evaluate, as the number of coefficients grows 
exponentially. Its main advantage is that it allows to fine-tune the number and characteris- 
tics of the paths, e.g. the maximum distance between end points. Fig. |6| shows the resulting 
errors when the error ellipsoid is restricted to a certain number of paths. It turns out that 
the relative error has a power law dependence on the number of paths, with an exponent of 
approximately —0.55 in this case. The contributions from different operators are relatively 
uniform except for the operators in the tails that contribute little or nothing, so the accuracy 
of the approximate Schur complement can be fine-tuned by varying the number of paths in 
the operator. 

The relative error measures how well the approximate Schur complement reproduces a 
delta function when applied to a Green's function. For actual applications, the more relevant 
quantity is the relative inversion error that looks at how good the inverse of the approximate 
Schur complement reproduces the Green's function of the original matrix. This is quantified 
by 

€ v = j:\fV -S-\U^)a^\ 2 (66) 

n 

with the relative error taken relative to the norm of f\ . This quantity was evaluated 
numerically by computing Green's function of the approximate Schur complement on the 
coarse lattice and comparing the results to the Green's functions of the original Wilson-Dirac 
matrix. Fig. [7] shows the relative error as a function of the approximation order both for 
the diagonal and the full optimization. It also decreases exponentially but is about twice 
the magnitude of the original error. 
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E. Spectrum 



The low eigenvalues of the original matrix M have a direct relation to those of the Schur 
complement Sm- If / is an eigenvector of M with eigenvalue A, 



Mf = Xf , (67) 



this becomes using (|16|) 



Sfr = A/i - M 21 M 22 1 Xf 2 
M 22 f 2 = Xf 2 - M 21 f\ . (68) 

Eliminating f 2 in the first line yields a modified eigenvalue equation for S: 

Sfr = A (l + M 12 M 2 - 2 \M 22 - \y x M 2 ^\ h . (69) 



If the eigenvalues of M 22 are large, as argued in sec. [Ill D| , the additional term on the 
right-hand side multiplying the eigenvalue is approximately 

1 + M 12 M 2 - 2 2 M 21 = 1 + Q+Q . (70) 

Assuming further that M 22 is short-ranged, the second term can be approximated by a 
constant a when acting on a low eigenvector, and (BUI) becomes 



Sh^\{l + a)h . (71) 

This rescaling of the eigenvalues accounts for the different lattice constant of the coarse 
grid. The appearance of the interpolation kernel Q in ( |70|) shows that it accounts for 
the parts of the original eigenvector that were on the fine grid and therefore dropped. In a 
geometric multigrid method, this is usually done by the coarsening prescription which moves 
information from fine to coarse degrees of freedom; in algebraic multigrid, it results directly 
from the LU decomposition. 



Fig. [10] shows the lowest positive eigenvalues of the Dirac matrix and different forms of the 
Schur complement for a sample U(l) gauge configuration at /3 = 3.0. With a rescaling factor 
of approximately 1 + a ~ 4, the exact Schur complement Sm reproduces the low eigenvalues 
of M well, and so does the order-2 optimized operator. The truncated von Neumann series, 
however, deviates significantly from the true results for the lowest eigenvalues, which are 
expected to be very important for the convergence properties when used as a preconditioner. 
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Relative error vs. number of paths 




# of paths in operator 



FIG. 8: The relative error as a function of the number of paths in the approximate Schur comple- 
ment for different orders of the path set. The paths have been reordered so that the relative error 
decays fastests. 

Rel. inversion error vs. max distance 
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FIG. 9: The relative error from the full optimization as a function of the maximum distance of 
path end points in the path set, shown for different maximum lengths of the path. 
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FIG. 10: Lowest positive eigenvalues of the effective coarse-grid operator S. The lower curve (filled 
circles) shows eigenvalues of the full Dirac matrix on a 16 x 16 lattice with a sample U(l) gauge field 
at /3 = 3.0. The three upper curves are the corresponding eigenvalues on the coarse grid, calculated 
from the exact Schur complement (full squares), the order-2 optimized operator (open squares), 
and the order-2 von Neumann series (open diamonds). The inset enlarges low eigenvalues. 



V. APPLICATIONS 



While a full discussion of the performance of multigrid operators in different algorithms 
is outside of the scope of this paper, we give here two examples of how they can be applied 
in numerical algorithms. Preconditioning uses the block LU decomposition provided by the 
multigrid decomposition to improve the condition number and thus convergence in iterative 
algorithms. Multigrid relaxation is the classical multigrid algorithm for solving the system 



(14]) and allows a comparison between a two-grid iteration and the original fine-grid iteration. 
In all calculations, a 16 x 16 lattice was used to allow the calculation of complete spectra. 
Note that only two-grid algorithms are considered; for realistic problems, the procedure 
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should be repeated on several grid levels. 



A. Preconditioning 

In preconditioning, the linear system 

Mf = a (72) 

is replaced by an equivalent system 

M~ x Mf = a' = M~ x a (73) 

using a preconditioning matrix M that is an easily invertible approximation to M. If M 
is sufficiently close to M, the spectrum of the preconditioned matrix M~ l M is contracted 
towards unity and its condition number reduced. For multigrid preconditioning, one uses 
an approximate block LU decomposition of the form ( |19D as the preconditioner. The pre- 
conditioned matrix is then 

/ 



M- l M 



1 \ I S- 1 1 Q + \ I M n M 12 
\ Q 1 ) { f - 1 ) [ 1 ) [ M 21 M 22 



(74) 



If Q and T are chosen to be exact, Q = — M 22 1 M 2 i and T = M 22 , this reduces to 



M~ X M 



I Q— 1 C 

d dm 



V 




(75) 



QS-^S-'Sm- 

The properties of the preconditioner are therefore determined by 1 — S^^-Sm, which is the 
quantity that was optimized in the optimization process fl4"T| ) above. 

Fig. [D] shows sample spectra of 1 — M~ 1 M with the Schur complement approximated to 
order 1, 2, and 3. If an exact Schur complement was used, the spectrum would reduce to a 
single point at the origin; the radius of the disk on which the eigenvalues lie in the complex 
plane is a measure of the quality of the preconditioner. The figure shows the result for two 
different choices of the interpolation matrix Q: In the bottom row, the exact interpolation 
Q = —M 22 M 2 i was used. Since Q is, as opposed to 5*, never inverted in the process, this is 
numerically legitimate and no more complex than the application of the approximate S -1 . 
The top row shows the result when an approximate interpolation matrix Q was used that 
was calculated in the same way as S in an optimization process. It shows that only if a 
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FIG. 11: Spectra of the multigrid preconditioned Dirac matrix 1 — M _1 M for a sample gauge field 
configuration. Optimized approximations of order 1, 2, and 3 where used from left to right for the 
Schur complement. The top row shows the spectra obtained using an approximate interpolation 
operator, the bottom row the spectra using the exact interpolation operator. Preconditioning with 
an exact Schur complement and interpolation operator would result in a single point at the origin. 
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FIG. 12: The preconditioned spectrum of fig. 11, using the von Neumann series as approximate 
Schur complement. 
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high order of approximation for S was chosen, it makes sense to use the exact interpolation 
operator, otherwise it just leads to a larger concentration of eigenvalues without reducing 
the radius of the disk. 

For comparison, fig. |T2| shows the spectrum of a matrix that was preconditioned using von 
Neumann series (and the exact interpolation operator) to the same order as before. While 
most modes are quite well reduced, there remain even at order 3 a few eigenvalues far away 
from the origin. These probably correspond to low eigenvalues of the original matrix that 



are not properly approximated by the von Neumann series, as seen in fig. 10 above. 



B. Multigrid relaxation 

Relaxation schemes for solving eq. (|nj) make use of an iterative process with an update 
step that is derived from the residual 

r M = a - M/ (n) , (76) 

where is the current approximate solution. The true error, i.e. the amount by which 
has to be updated, can be calculated from the residual by 

e (n) = f _ f(n) = M -l r W _ (77) 

Given an approximation M _1 to M _1 , an approximate update step reads 

f(n+l) = f(n) + Tr (n) = ^ _ jtf-ljtf^ (n) + Tfl _ (73) 

This iteration always has a fixed point at the true solution /, and the quality of the approx- 
imation M~ l to M _1 only determines the convergence characteristics of the process. 

Multigrid relaxation is based on using the Schur complement along with restriction and 
interpolation to construct an approximation to M _1 on the coarse lattice. First, the residual 
is restricted to the coarse lattice using the restriction 

fW= ei + Q + e 2 . (79) 

A coarse lattice approximation e^ 1 ' to the true error is then found by applying the inverse 
of the approximate Schur complement: 

g(n) = £-l f (n) _ (g ) 
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FIG. 13: Spectra of the convergence matrix for multigrid relaxation in a sample gauge field. Opti- 
mized approximations of order 1, 2, and 3 where used from left to right for the Schur complement. 
The top row shows the convergence of the multigrid cycle if a full inversion is performed on the 
coarse grid, and one relaxation step on the fine grid. The bottom row shows the same spectrum 
if only a single relaxation step is performed on each grid. The convergence rate is given by the 
radius of the smallest disk containing eigenvalues. 
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Finally, this is interpolated back to the fine lattice: 




(a) 




(81) 



and used to update the approximation. This update step acts on the coarse degrees of 
freedom only and must be followed by a compatible relaxation step on the fine lattice. In 
matrix form, the approximate inverse of M used here reads 



It differs from the true LU decomposition (p0|) in that the lower right entry of the matrix 
in the middle is zero instead of M^ 1 . If a perfect interpolation and restriction is chosen, 
Q = — M 2 ^ 1 M2i, the iteration matrix for the coarse-grid step reads 



Its convergence properties on the coarse lattice are thus given again by 1 — S~ x S as in the 
case of preconditioning; on the fine lattice, it performs no relaxation at all, and must be 
followed by an ordinary relaxation step using M 22 . 

Fig. [13] shows the spectra of the convergence matrix for multigrid relaxation. In the top 
row, each step consists of performing a complete relaxation on the coarse lattice, and one 
relaxation step on the fine lattice. This leads to a concentration of eigenvalues around zero, 
showing modes that are nearly completely reduced on the coarse grid. The remaining modes 
must be reduced on the fine grid and determine the residual convergence rate. In an actual 
application, one would not completely reduce the modes on the coarse grid, but perform 
coarse- and fine-grid relaxation alternately. The resulting convergence spectrum is shown in 
the bottom row; it has a similar convergence rate as the complete reduction, but requires 
much less operations. Note that increasing the order of the approximation from 2 to 3 does 
not improve the convergence, but this might be related to the fact that only one fine-grid 
relaxation step was used here. 

Finally, we compare the number of operations required to calculate a sample Green's 
function in U(l) lattice gauge theory on a 16 x 16 lattice. It must be cautioned that this is 
only an example calculation using two-grid relaxation, not a full multigrid algorithm. The 
method offers a multitude of choices and tunable parameters: the order and type of the 




(82) 




(83) 
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Two-grid cycle 1 :1 
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FIG. 14: Sample convergence of a two-grid relaxation method compared to the same method 
applied to the base grid, for a U(l) Dirac operator on a 16 x 16-lattice. The horizontal axis gives 
the number of floating-point multiplications, the vertical axis the true error of the algorithm. For 
the two-grid algorithm, one sweep on the coarse lattice and one sweep on the fine lattice were 
performed alternately. 



approximation for the LU decomposition, the iterative algorithm (used here is for simplicity 
Jacobi relaxation; in actual applications one would probable use Gauss-Seidel), the multigrid 
cycle, and the convergence criteria on the two grids. The size and dimensionality of the 
lattice and the gauge group, as well as the amount of disorder in the gauge field, have been 
seen to decisively influence results in other algorithms. Also, the two-grid algorithm must 
be extended to a true multigrid algorithm by recursively applying the method. And finally 
relaxation might be less performant than a preconditioned Krylov solver. The ultimate 
measure of performance will be the actual computer time used in the algorithms, but this 
will depend on the implementation and on the architecture of the target machine. For 
example, the computation of the approximate Schur complement can either be performed 
beforehand and stored in memory, or on the fly in the matrix-vector multiplication of the 
Schur complement, depending on the amount of memory available. Here a lower order 
of approximation, resulting in a smaller Schur complement, can be favorable though its 
convergence rate per relaxation step might be slower. 
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With these precautions, an example run for calculating a sample Green's function is 
shown in Fig. [TJ]. As a measure of computer time used, the approximate number of floating- 
point multiplications is given. The convergence is measured as the error norm against the 
true result. 

VI. CONCLUSIONS 

It was demonstrated how the algebraic multigrid method can be applied to disordered 
linear operators. Other than in previous multigrid approaches, a thinned lattice for the 
coarse degrees of freedom is used and averaging over block spins is avoided. The interpola- 
tion operator is chosen to obtain a block LU decomposition, and the resulting coarse-grid 
dynamics is completely described by the Schur complement of the original matrix. Gauge 
covariance is automatically taken into account by this procedure. 

The Schur complement can be approximated to arbitrary precision in a numerical opti- 
mization process by expanding it in a linear basis constructed from connected paths on the 
original lattice, forming a generalized stencil for disordered operators. Similar to renormal- 
ization group techniques, the numerical procedure for obtaining the expansion coefficients 
uses the coarse-grid projection of an ensemble of systems to obtain the coefficients of the 
approximate coarse-grid operator. In this way, the information gathered in the optimization 
process can be used to speed up the actual calculations as the optimized effective coarse-grid 
operator "learns" the dynamics of the system. 

The resulting effective coarse-grid operator is constructed from the generalized stencil as 
a weighted sum over paths on the original lattice. It forms a denser, but smaller matrix with 
next-neighbor and higher interactions and can be used to improve performance in numerical 
algorithms by preconditioning or similar methods, hopefully not only for matrix inversion, 
but also for other problems such as eigensystem analysis and rational matrix functions. A 
variety of tunable parameters make it possible to adapt the procedure to different algorithms 
and architectures. Whether the method will actually improve real-life algorithms, remains 
of course to be seen until realistic systems and efficient implementations are investigated. 
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